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Abstract. The upper critical field in MgB2 is investigated in the framework of the two-gap Ginzburg- 
Landau theory. A variational solution of linearized Ginzburg-Landau equations agrees well with the Landau 
level expansion and demonstrates that spatial distributions of the gap functions are different in the two 
bands and change with temperature. The temperature variation of the ratio of two gaps is responsible 
for the upward temperature dependence of in-plane H C 2 as well as for the deviation of its out-of-plane 
behavior from the standard angular dependence. The hexagonal in-plane modulations of H C 2 can change 
sign with decreasing temperature. 

PACS. 74.70. Ad Metals; alloys and binary compounds (including A15, MgB2, etc.) - 74.20.De Phe- 
nomenological theories (two-fluid, Ginzburg-Landau etc.) - 74. 25. Op Mixed states, critical fields, and 
surface sheaths 



1 Introduction 

Multigap superconductivity PE] has been discussed in 
the late 1950's for materials with a varying strength of 
electron-phonon interactions between different pieces of 
the Fermi surface. After the discovery of superconductiv- 
ity in MgB2 [3] in 2001, an impressive collection of exper- 
imental and theoretical works 0] has established that this 
compound is the first unambiguous example of a multi- 
gap superconductor. In MgB 2 the charge carriers are dis- 
tributed between two sets of bands: the cr-bands with 
quasi-2D cylindrical Fermi sheets and the 7r-bands with 
3D sheets forming a tubular network. The electron-phonon 
coupling is stronger in the cr-bands than in the 7r-bands, 
and gives rise to an s-wave phonon-mediated supercon- 
ductivity with two gaps A± ~ 7 meV and A 2 ~ 2.5 meV. 
Since the two sets have different characteristics (interac- 
tion with phonons, geometry of the Fermi sheets, impurity 
dependence etc.), an interplay between them results in de- 
viations from the standard BCS theory. The most striking 
consequences of the two gaps are the unusual anisotropic 
features of MgB 2 under magnetic field, for example, in- 
equality between the penetration depth and the upper 
critical field anisotropies, and their variations with tem- 
perature ppir7ll51PlTnillllll^lll3lll4lll5lll6| . and the 30°- 
reorientation of the flux line lattice with increasing mag- 
netic field applied along the c-axis |17II18| . 

The two-gap Ginzburg-Landau (GL) theory for MgB2 
developed in Ref. JS] (see also the preceding works ^3 
I20| ) is the exact limit of the microscopic theory in the 
vicinity of the transition temperature. It can thus account 
for most of the observed properties in a clear and coher- 



ent way near T c , while its simplicity compared to earlier 
studies is useful to understand the physics in this mate- 
rial. In the present paper we extend our previous analysis 
of the two-band effects ^H] on angular and temperature 
dependence of the upper critical field H c2 . We minimize 
the GL functional using a variational procedure, which 
highlights separate spatial anisotropies of the gap in each 
band. This is an improvement compared to the earlier so- 
lutions where only one common distortion for both gaps 
is considered |11U14| . This method is compared to a solu- 
tion based on the Landau level expansion. We then esti- 
mate the temperature range of the GL regime. The present 
study covers the out-of-planc H c2 anisotropy. By going be- 
yond the ellipsoid Fermi sheet approximation of Refs. 
114) . we also calculate in-plane modulation of the upper 
critical field arising from the hexagonal crystal symmetry. 

For a clean two-band BCS superconductor with two 
gaps A\ and A2, the GL functional ^2] has the form 



F, 



GL 



dx 



ail^il 2 + a 2 \A 2 \ 2 - ~f(AlA 2 + A 2 Ai) 



K li \n i A 1 \ 2 + K 2l \n t A 2 f + i/3i|Zii| 4 + i/3 2 |zi 2 
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c 
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where repeating index i implies a sum, <Pq is the quan- 
tum flux, A the potential vector, gi and g 2 the intraband 
pairing coefficients (n=l,2 for the cr, 7r-band), 53 the inter- 
band pairing coefficient, G — gig 2 — gf, N n the density 
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of states at the Fermi level in the band n, u>d the Debye 



frequency, C the Euler constant, and 



the square of 



the Fermi velocity i-component averaged on the sheet n. 
It is then convenient to write u\ — —a\t with a\ = N%, 
t = ln(Ti/T) and T x = (2u D e c /Tr)e- S2 / GNl for the first 
active band and oti — ct2o — &2t with ai = N2 for the 
passive band. In MgB2 the active and passive bands cor- 
respond to the er and n bands, respectively. 

The crystal structure of MgB2 is uniaxial, so the gra- 
dient coefficients are the same for all directions in the 
basal plane, and K na = K no = K n . LDA calculations [2T] 



yield for the highly anisotropic er-band (v 



Flab/ 

(v'pic) — 0.05, while for the 7r-band (v F2ab ) — 1.51 and 
( v F2c) = 2.96, all numbers are in units of 10 15 cm 2 /s 2 . 
With the provided ratio N 2 /Ni = 1.5, the in-plane gradi- 
ent constants for the two bands are practically the same 
K2/K1 w 1.06, whereas the c-axis constants differ by al- 
most two orders of magnitude K2 C /Ki c w 90. A crude 
estimate for H C 2 at zero temperature by ^XiK 1^ ~ 
4T is substantially smaller than the experimental value 
-ff?|(0K) as 18T, which suggests the gradient constants 
based on LDA data are over-estimated by a factor of four. 
Such a discrepancy is due to a significant renormalization 
of effective masses by the electron-phonon coupling. The 
electron-phonon coupling leads to effective masses twice 
larger than the LDA prediction in the cr-band, whereas 
they are only slightly renormalized in the 7r-band f22] . The 
reduction of gradient term coefficients is given by squares 
of the mass renormalization factors. 

The interband impurity scattering in MgB2 is excep- 
tionally small due to its particular electronic structure, 
even in low quality samples |23| . The clean limit two-gap 
GL theory described above is straightforwardly extended 
to include the effect of s-wave intraband scattering by non- 
magnetic impurities |24| : the GL functional keeps the same 
form wherein the expression for K n i has to be replaced by 
K ni = nN n (v Fm )A{T n )/8T c with 



2.13 and 



A(r n ) 



7T 2 ^ (2m + l) 2 ((2m + 1)27tt„T c + 1) 



1 



(2) 



where t„ is the transport collision time in the band n. The 
intraband anisotropy is then the same as in the clean limit, 
while the renormalization factor A(r n ) can vary between 
the two bands due to different sensitivity to impurities. 
The resulting GL equations are naturally found as the 
limit of Usadel equations near T c . |lllll5j . 



2 Upward curvature of H$(T) 

In this section the z-axis is fixed along the crystal c-axis 
and the y-axis is taken parallel to the magnetic field ap- 
plied in the afe-plane. The vector potential is chosen in the 
Landau gauge as A = (Hz, 0,0). The coupled linearized 
GL equations for solutions homogeneous along the field 
direction are 



K n h 



V 



K nc dl) A n - 7 A 







(3) 
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Temperature, K 



Fig. 1. The upper critical field in MgE^: experimental data 
from Lyard el al. (symbols) and the GL computations with 
parameters 020 = 0.65ai and 7 = 0.4ai (solid lines). The inset 
shows the anisotropy ratio = H^/H^. 



for n — 1,2, n' = 2,1, with the reduced magnetic field 
h = 2nH/(p . Since K lc /Ki ^ K 2c /K2, an analytic solu- 
tion can not be obtained by rescaling distances as in the 
single-gap case. We, therefore, search for an approximate 
solution of the form 



4i 

±2 



(4) 



where the Landau level wave functions are defined by 



1 



1 
71 



{-d z + z) 



-z 2 /2 



(5) 



Different coherence lengths for each band are allowed with 
the parameterization £ 2 = jj, n /h where /i„ quantifies the 
distortion of the spatial distribution of the n-th compo- 
nent (in the single-gap case, fj, is the stretching factor of 
the flux line lattice at the upper critical field and is inde- 
pendent from temperature). The following quadratic form 
in the GL functional is then found: 



F 2 = (ai +hK 1 )\c\ 2 
-j(c*d + d*c) 



(a 2 + hK 2 ) 



(G) 



with K n - ±{K nt i n + K nc /fi n ) and 7 2 = 7 2 ^gf ■ At 



the transition field, the determinant in Eq. © vanishes. 
This condition leads to 



ai 
2Ki 




a 2 

2Kn 



o 2 



K t K 2 
(7) 

In order to find the (nucleation) upper critical field, h is 
maximized h c2 = iaax.^ ljfl2 h(fii,fi 2 ). 

Within the above variational scheme, the analytic ex- 
pressions for the in-plane transition field are possible in 
two temperature regimes. Near T c , vanishing h implies 
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Fig. 2. In-plane H C 2 calculated with the Landau level expan- 
sion to the order iV max for the same parameters as in Fig. 
The inset displays the gap ratio p = A2/A1 found with the 
highest expansion order. 



that the superconducting gaps have the same variation 
length in each band. The condition fi± = \i% yields 



n c2 



a\{t - t c 



yJ{K x + p 2 K 2 ){(K Xz + p 2 K 2z ) 



(8) 



with the gap ratio p = \d/c\ « 1/01/02 and t c = ln(Ti/T c ). 
bmce Kic ~ 0.0Lff 2c in MgB 2 , whereas p 2 « 0.1, we can 
simplify the above expression to 



h ab 

n c2 



a\{t - t c 



py/K 1 K 2 



(9) 



In the second temperature regime for T < T±, the first 
active band is dominant and 



lab 
n c2 



ai t 



(10) 



lc 



The line if"|(T) exhibits, therefore, a marked upturn cur- 
vature between the two regimes, in contrast to H£ 2 (T). 
The two upper critical fields are plotted in Fig. ^ In order 
to fit the experimental data, we have renormalized all gra- 
dient constants obtained from the LDA data by a factor 
of five. The corresponding mass enhancement y/b w 2.2 
roughly agrees with the electron-phonon renormalization 
factor 22] . For simplicity, the same value has been applied 
for both bands. 

In order to verify an accuracy of the variational method, 
we alternatively proceed by expanding the gap functions 
in terms of the the Landau levels: A n = ^2 p c n , P 4>n,p where 

= & 1/2 *p(z/t n ) and £ = y/K nc jK n /h. For the 
upper critical field this expansion is restricted to the even 
order levels. The quadratic part of the GL functional has 
the following matrix clement in this base: 



with n,n' S {1,2}, and p, q > 0. The upper critical field 
h c2 is then approximated by the largest root of the sub- 
matrix determinant corresponding to the desired expan- 
sion up to the order iV max . 

Although the zeroth order approximation significantly 
deviates near T c (see Fig.|2J, the procedure is rapidly con- 
verging with increasing the expansion order, even in the 
case of a great disparity between the two bands (e.g., 
> 100 or < 0.01). The expansion to the order 
A^ m a X > 12 yields the upper critical field curve in ex- 
cellent agreement with the variational solution (the two 
curves are indistinguishable on the scale of Fig. |2J. 

Fig.[2]displays the behavior of the parameters p n defin- 
ing the effective anisotropy of the variation lengths £ n 
in the plane perpendicular to the magnetic field, i. e. 
fi>n = Cnc/^nab for the magnetic field applied in the basal 
plane. This confirms the above analytic predictions: the 
order parameter varies on different length scales for each 
band, and \i n can change with temperature contrary to 
the single-gap case. At T c , the two parameters have the 
same value y/(K lz + p 2 K 2z )/(K 1 + p 2 K 2 ) = 0.59 with 
p = 0.44, while p, x « y/K^/Wi = 0.15 below T x = 29 K. 
We should stress that periodic vortex structures for the 
two gaps have the same lattice parameters for arbitrary 
ratio of pi/p 2 . However, spatial distributions of |Z\i(r)| 2 
and |Z\ 2 (r)| 2 become quite different at low temperatures 
once p\ -C p 2 . Such a behavior is demonstrated on the 
top panel of Fig. The different spatial distributions of 
the two gaps can be probed by scanning tunneling mi- 
croscopy. Also, magnetic field generated by superconduct- 
ing currents h s (r) ~ (|Z\i(r)| 2 + |zA 2 (r)| 2 ) should deviate 
significantly for a distribution expected for an anisotropic 
single-gap superconductor. Muon spin relaxation measure- 
ments can in principle verify such a behavior. 

We shall now estimate the temperature range of the 
GL regime from the above computations. The gradient ex- 
pansion is valid as long as \K n idfA n \ < A n for all n and i. 
This condition is approximately replaced with K n ijQ li < 
1. The most restrictive case is for K 2c j^ c — K 2c h^/ p 2l 
which becomes ~ 1 below ~ 30 K, well beyond a nar- 
row temperature regime suggested for the GL theory by 
Golubov and Koshelev fTJj. The discrepancy is partially 
terminological, since in Ref. JI] the GL approximation al- 
ways corresponds to an effective (anisotropic) single-gap 
GL theory, which is correct only when the ratio of the two 
gaps is constant. As we have demonstrated above, the full 
two-gap GL theory is valid in a much wider temperature 
range and describes adequately temperature variation of 
A 2 /A\ (Fig.EJ and of the two coherence lengths (Fig. EI) - 



3 Angular dependence of out-of-plane H, 



c2 



M F2 

lvl 2p+n,2q+ri 



"7 dx 4>* nap {x)4>r, 



n.n' Or 

,2q(x) (1 - 5 n ,n') 



(11) 



Let us now discuss the out-of-plane behavior of the upper 
critical field. In the single-gap anisotropic GL theory, when 
H is tilted from the c-axis by an angle 9, the upper critical 



4 



V. H. Dao, M. E. Zhitomirsky: Anisotropy of the upper critical field in MgB2 



IAJ at 30K 




|A J at 30K 



|A, | at ~35K 



CD 

H — ' 

CD 

E 

CO 
i— 

CO 

Q- 

"cG 

c 
o 

-i — ' 

CO 
CO 

> 



0.6 - 



0.0 




21 28 

Temperature, K 



35 



Fig. 3. Lower panel: temperature dependence of variational 
parameters for in-plane magnetic field and the same set of GL 
parameters as in Fig. Upper panel: absolute values of the 
two gaps in the vicinity of H C 2 near to and away from the 
transition temperature. 



field has an elliptic (effective mass) angular dependence 

H c c2 (T) 



H SAGh (6,T) = 



cos 2 (6) +sin 2 (6>) 7 ^ 2 



(12) 



where jh = H^/H^ is a temperature independent con- 
stant \J Kc/Kab. Experimental measurements in MgB2 
have shown that not only "/h changes with temperature 
(Fig. ^) but deviations from the elliptic angular depen- 
dence (I12|) grow with decreasing temperature [3]|HK|. Such 
a behavior has been reproduced within quasi-classical Us- 
adel equations ^TJ- The methods we have employed for 
arc still valid to find H C 2(8): one needs only to re- 
place K nc by an angular dependent K n (6) = cos 2 (9)K n + 
sin (6)K nc in the previous formula. Expression (0) for h 
shows that the deviation grows with the disparity be- 
tween the K n {6), so it increases when departing from 
T c . The deviations can be quantified by 5A(9) = 1 — 
(H c2 (9) I H S agl(0)) 2 ■ Fig. ^ displays the maximum devi- 
ation SA max = ma,xgSA(9). The dashed line is obtained 
from the two-gap GL theory with the parameters used 
above to fit the f/^-data by Lyard et al.]§\ in Fig.f^ The 
calculation qualitatively reproduces experimental data from 
Rydh et al. (Hj: 5A m&x increases with decreasing temper- 
ature and then saturates. But a quantitative discrepancy 
appears below 0.9T C and becomes important at lower tem- 
perature. This deviation can be partially explained by the 
fact that experimental results are strongly sample depen- 
dent. At the present, the origin of the discrepancy remains 



< 

CO 




0.9 1.0 



Fig. 4. Maximum deviation from the single-gap GL scaling 
law : squares are experimental data [§] and lines are obtained 
from the two-gap GL theory. 



an opened question. The full line is obtained with a mod- 
ified interband coupling 7* = 0.3ai and K* 2 = if2c/3 
corresponding to a smaller anisotropy in the 7r-band. 



4 In-plane modulation of H, 



c2 



In a hexagonal crystal, the transition magnetic field should 
exhibit a six-fold modulation when rotated about the c- 
axis j2S]. The crystal field effect on superconductivity can 
be incorporated to the GL theory by including higher or- 
der (non-local) gradient terms |2f>| . Symmetry arguments 
suggest that coupling between the superconducting order 
parameter and the hexagonal crystal lattice appears at the 
sixth-order gradient terms. For a two-gap superconductor 
like MgB2, the additional sixth-order part of the free en- 
ergy is a sum of separate contributions from each band: 
Fql — ► Fgl + ^6,1 + ^6,2- The correction derived from 
the BCS theory [181127) is (omitting the index n = 1, 2 for 
brevity) 



Fa 



C(7)N 



1 



1 

2^ 



32tt 6 T c 6 

x {n l n J n k A)*(iiiii m iinA) 



V Fi V Fj V Fk V F lVFmVFn) 



(13) 



Setting the z-axis perpendicular to the basal plane, the 
above terms can be split into isotropic in-plane part 



F™° = K l 6 A" 



n 2 + n 2 

x 1 y 



(14) 



withal = ^(1-^)(( W 
in-plane contribution 



( v Fy)) ' an d anisotropic 



FS 



1 



KSA* 



(n x + m y f + (n x ~ in y f A (15) 
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5. Influence of the isotropic correction F$ with 
, («| 2 » = i/(4.6, f .6) (in units of 10 46 (cm/s) 6 ) for dif- 

magnitudes v: relative shift of the upper critical field 
Hc 2 — H° 2 )/H° 2 and corresponding variational parame- 



-k 6 , 2 h 3 )\d\ 2 



j(c*d + d*c) 



(17) 



with Kq iTI 

this results in the new matrix element 



= -^ii"6,nM„- While in the expansion method, 



M. 



F 2 +F 6 
2p+n^2q+n' 



— mF 2 



(18) 

with Mpfg = -| (^ p |(a t + a) 6 \& q ) where a is the annihila- 
tion operator of Landau levels. 

In the weakly anisotropic regime Fq « F 2l we expect 



H c2 (<l>) « H° c2 (1 + rf + r? a cos(6</>)) 



(19) 



The isotropic parts yield a ^-independent shift of H c2 (and 
ensure K§ n > for the numerical solution converging) 
while the anisotropic parts are responsible for the six-fold 
modulation of the correction, rf can change sign when the 
temperature varies because the anisotropics in each band 
are opposite. Fig.Eldisplays the corrections brought by the 
isotropic parts of jpg. The deviations become important 
below 30K as expected out of the estimated GL regime, 
which implies the necessity to retain higher order terms 
in the gradient expansion of the GL functional. 

The extra h 3 terms prevent from deriving an analyt- 
ical expression for the magnetic field correction 5h c2 — 
h C 2{4>) — /i C 2(y/12). We can however partially estimate the 
latter. Let us name the quantities related to the quadratic 
form (F 2 + Fg s °) with the superscript "o", and the ones 
for (F 2 + F^ so + F 6 an ) without it. Within the variational 
method, we then find with a perturbation expansion 



with Kt = (1 - «t,« x ) - (v Fy )). This expres- 

sion of Fq 11 assumes that the x- and the y-axes are parallel 
to the reflection lines in the a6-plane. With the x-axis par- 
allel to the 6-direction, tight-biding calculations ^2] yield 



J Fxl 



4.608, (v 



the 7r-band, 



Fyl 
= 1.514, 



4.601 for the c-band, while for 



J Fy> 



1.776 in units of 10 



Mi 



(cm/s) 6 . The different sign of the hexagonal harmonics of 
the Fermi velocities in the two bands is responsible for 
a unique 30-degree orientational transition of the vortex 
lattice in MgB 2 ._18. No theory can describe at present the 
elcctron-phonon effect on the hexagonal modulation of the 
Fermi surface. We use, therefore, the raw LDA values for 
all gradient coefficients in the consideration below. If we 
rotate now the orthogonal axes so that the y-axis is paral- 
lel to the magnetic field H when the latter forms an angle 
<fi with the a-axis, the terms in Fq change in a simple way: 
is preserved while F§ n turns into 



771S0 

^6 



pan 
^6 



-KSA* 



\n x + m y f + e- t6 *(n x ~in y 



A 

(16) 

Since FI y = 0, the extra term can be written as Fq = 

Kq n cos(60). For 



6,n 



Y,Ke,nA^II%A n with Ke, n = K, 
the variational approximation, the new functional yields 
the quadratic form 

F 2 = {at + K x h + K 6A h 3 )\c\ 2 + (a 2 + K 2 h 



Sh 



c2 



1 dh 



c2 



(Z\°|F 6 an |zA° 



h° dT (A°\d(F 2 +F™)/dT\A° 



'<:2 

15 



(20) 



J — -r^n r0 ■ - 



dT 



a>i + a2P 



where p 2 = \c° /d°\ 2 . The expansion method provides in a 
similar way 



Sh c2 T dh° c2 (l + p 2 )(A°\F£ n \A° 



>c2 



K2 dT 



ai + a 2 p A 



(21) 



but (Z\°|.F 6 an |Z\ 



has a more complicated expression. 
In Fig. El we have plotted the relative modulation 
amplitude r, a = (/i c2 (0) - h c2 (-ir/6))/(h e2 (0) + K 2 {-k/&)) 



with the hexagonal anisotropy (v F y 
(v%„) = (vlY - {v%) a where ((v F1 ) 



J F2l 



v%) a and 
(4.6,1.6) 

46 



while ((v F1 ) a , (v 6 F2 ) a ) = (0.3^,-0.2^) (in units of 10 
(cm/s) 6 ). Ab initio calculations provides (v F1 ) a <C (v F2 ) a 
for MgB 2 which corresponds around to the couple (y\, v 2 ) = 
(0; 1) in Fig. |SJ Due to the LDA results uncertainty and 
also to illustrate the interplay between the two bands, the 
plots for other values of (vi,v 2 ) are displayed. Note the 
results at low temperature should be taken with caution 
since they are obtained out of the GL regime. When the 
hexagonal anisotropies of each band are of the same or- 
der, rf sign can change with temperature. But this mod- 
ulation is too small to be detected experimentally in the 
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Fig. 6. Relative modulation amplitude for different pairs 
(yi,va)'- solid lines are obtained with the variational method 
and symbols with estimate 



GL regime, which agrees with measurements reported by 
Shi et al. [Hj- Estimation (|20() gives three reasons for this. 
First, rj a grows as h 2 2 contrary to the four- fold symme- 
try crystal case where the increase is linear. Then the 
anisotropics of the two bands oppose each other. And fi- 
nally, even though (v® Fl ) a would be too small to compete 
with (vp 2 ) a i the contribution from the second band is re- 
duced by the rapidly decreasing factor p 2 and, below 30K, 
by (4- 



5 Conclusions 

Angular and temperature dependence of the upper critical 
field of MgB2 have been determined within the two-gap 
GL theory. We have used two different numerical meth- 
ods which are in excellent agreement with each other and 
yield an unconventional anisotropy of H C 2 observed in the 
superconductor MgB2. Such a behavior reflects the differ- 
ent Fermi sheet geometries and the varying importance 
of the small 7r-gap. The zeroth Landau levels employed 
in the variational approach are sufficient for accurate de- 
scription of the continuous transition at H C 2- Contrary to 
the single-gap case, spatial anisotropy of the gap func- 
tions in the plane perpendicular to the magnetic field 
changes with temperature and can be different for each 
band. This explains the deviation from the effective mass 
angular dependence (jf 2JI applicable to ordinary supercon- 
ductors. Existence of two different characteristic lengths 
should also affect the vortex core shape, ^B] especially 
when an applied field is perpendicular to the c-axis. The 
gap functions have an effective single-component behavior 
only in a temperature region near T c significantly narrower 
than the range for the validity of the two-gap GL theory 
~ (T c - T)/T c ~ 1/7. At last, the hexagonal afe-plane 
modulation of H C 2 arising from the crystal symmetry can 
result in a change of the sign of the hexagonal harmonics 
of H C 2(9) when the temperature is decreased. 
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